rm(list=ls())
library(data.table)
library(ggplot2)
library(polywog)
library(MASS)
library(doParallel)
library(directlabels)

load("ensemble.Rdata")
ensemble$sddistance <- ensemble$sd.distance/.12
ensemble <- ensemble[, sd.distance:=NULL]
ensemble <- ensemble[, diff:=parties-mean.num.gov.parties]
ensemble <- ensemble[,alpha5:=ifelse(alpha>.5, 1, 0)]


#########################
# FUNCTION FOR PLOTTING #
#########################
plot.vars <- function (data, varying, by="",
                         parties.val=5,
                         discount.val=.5,
                         alpha.val=.5,
                         sd.distance.val=1
) {
 if (varying=="parties"|by=="parties") {
     parties.val <- unique(data[,parties])
 }
 if (varying=="discount"|by=="discount") {
     discount.val <- unique(data[,discount])
  }
 if (varying=="sd.distance"|by=="sd.distance") {
     sd.distance.val <- unique(data[,sd.distance])
 }
 if (varying=="rule"|by=="rule") {
     rule.val <- unique(data[,rule])
 }
 if (varying=="alpha"|by=="alpha") {
     alpha.val <- unique(data[,alpha])
 }
 output <- data[round==76&
                   discount %in% discount.val
                 & parties %in% parties.val
                 & alpha %in% alpha.val
                 & sd.distance %in% sd.distance.val]
 output <- as.data.frame(output)
 plot.data <-data.frame(x=output[,varying],
                          y=output$mean,
                          lower=output$mean-(1.96*output$se),
                          upper=output$mean+(1.96*output$se),
                          variable=output$variable,
                          rule=output$rule)
 if (by!="") {
   plot.data$z <- output[,by]
 } 
     p <- ggplot(plot.data, aes(x=x, y=y, group=rule))
     #p <- p + geom_segment(aes(colour=variable, y=upper, yend=lower, xend=x))
     #p <- p + geom_point(aes(colour=variable))
     p <- p + geom_ribbon(fill="grey73", aes(ymin=lower, ymax=upper))
     p <- p + geom_line(aes(colour=as.factor(rule)))
     p <- p + theme_bw()
     p <- p + theme(panel.grid.major.x = element_blank())
     p <- p + theme(panel.grid.major.y = element_blank())
     p <- p + theme(panel.grid.minor.x = element_blank())
     p <- p + theme(panel.grid.minor.y = element_blank())
     p <- p + xlab(varying)
     p <- p + ylab("")
     p <- p + scale_colour_manual("Rules",
                                  breaks=c(2, 3, 4, 7, 8),
                                  values=c("red", "blue", "yellow", "green", "pink"),
                                  labels=c("Aggregator", "Sat Governator", "Sticker",
                                        "Hunter", "Hunt Governator"))
     ## p <- p + scale_colour_manual("Rules", labels=c("Aggregator", "Sat Governator", "Sticker",
     ##                                    "Hunter", "Hunt Governator"))
     if (by!="") {
         p <- p + facet_grid(variable~z, scales="free")
      } else {
         p <- p + facet_grid(variable~., scales="free")
      }
return(p)
} #end plot.vars

#This function data-mines using polywog but eliminates terms from the equation
polywog.select <- function(formula, data, degrees.cv, boot=5, select=.05, lambda=NULL) {
    cv <- cv.polywog(formula=as.formula(formula),
                     data=as.data.frame(data),
                     degrees.cv=as.numeric(degrees.cv),
                     .parallel=TRUE)
    fit <- cv$polywog.fit
    fit <- bootPolywog(fit, nboot = boot, .parallel=TRUE)  
    sig <- which((abs(summary(fit)$coefficients[,1]/summary(fit)$coefficients[,2])>1.96)==TRUE)
    rhs <- names(sig)
    rhs <- gsub("[.]", ")*I(", rhs) 
    if (rhs[1]=="(Intercept)") {
        rhs <- rhs[-1]
    } else {
        rhs[1] <- "-1"
    }
    if (length(rhs)>1) {
        rhs <- ifelse(rhs!="-1", paste("I(", rhs, ")", sep=""), "-1")
        rhs <- paste(rhs, collapse="+")
    } else if (length(rhs)==1) {
        rhs <- paste("I(", rhs, ")", sep="")
    } else if (length(rhs)==0){
      rhs <- "1"
    }  
    new.formula <- as.formula(paste(as.character(formula[2]), "~", rhs))
    return(lm(new.formula, data))
}

#This function makes predictions based on polywog.select results.
predict.polywog <- function(pmodel, parties=3, discount=0, sddistance=0, alpha=0, rule="") {
 input.data <- expand.grid(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  parties, discount, sddistance, alpha)
 names(input.data) <- c("mean.eccentricity", "mean.representation", "mean.gov.representation",
                        "mean.gov.representation2",
                        "mean.num.gov.parties", "mean.gov.seats", "mean.gov.cohesion",
                        "mean.oversized", "mean.majority", "mean.ENgov", 
                        "parties", "discount", "sddistance", "alpha")
 input.data <- as.data.frame(input.data)
 output.data <- as.data.frame(input.data)
 input.data$alpha1 <- ifelse(input.data$alpha==1, 1, 0)
 input.data$alpha0 <- ifelse(input.data$alpha==0, 1, 0)
 input.data$alpha5 <- ifelse(input.data$alpha==.5, 1, 0)
 input.data$discount1 <- ifelse(input.data$discount==1, 1, 0)
 input.data$discount0 <- ifelse(input.data$discount==0, 1, 0)
 input.data$sddistance0 <- ifelse(input.data$sddistance==0, 1, 0)                         
 input <- model.matrix(pmodel, data=input.data)
 betas <- mvrnorm(1000, coef(pmodel), vcov(pmodel))
 values <- input %*% t(betas)
 means <- apply(values, 1, mean)
 lower <- apply(values, 1, quantile, p=.025)
 upper <- apply(values, 1, quantile, p=.975)
 return(as.data.table(cbind(output.data, lower, means, upper,
                            type=as.character(pmodel$terms[[2]]),
                            rule=rule)))
}



#This function plots based on predict polywog.
plot.polywog <- function(x, varying, by="", plot=TRUE, color=FALSE, interest="eccentricity",
                         original=NULL) {
 interest.mean <- paste0("mean.", interest)
 interest.se <- paste0("se.", interest)
 x <- as.data.table(x)
 if (plot==TRUE) {
   x <- as.data.frame(x)
     if (by!="") {
       plot.data <- data.frame(x=x[,varying],
                               y=x$means,
                               lower=x$lower,
                               upper=x$upper,
                               z=x[,by],
                               original=0)
#      plot.data <- as.data.frame(cbind(x,y, y.se, z))
      } else {    
       plot.data <- data.frame(x=x[,varying],
                               y=x$means,
                               lower=x$lower,
                               upper=x$upper,
                               original=0)
   }
   if (is.null(original)==FALSE) {
       #original <- as.data.frame(original)
      if (varying=="parties") {
          parties.val <- unique(c(original[,parties], x[,"parties"]))
      } else {
          parties.val <- unique(x[,"parties"])
      }
      if (varying=="discount") {
          discount.val <- unique(c(original[,discount], x[,"discount"]))
      } else {
          discount.val <- unique(x[,"discount"])
      }
      ## if (varying==memory) {
      ##     memory.val <- unique(c(original[,memory], x[,memory]))
      ## } else {
      ##     memory.val <- unique(x[,memory])
      ## }
      ## if (varying==evol) {
      ##     evol.val <- unique(c(original[,evol], x[,evol]))
      ## } else {
      ##     evol.val <- unique(x[,evol])
      ## }
      if (varying=="sddistance") {
          sddistance.val <- unique(c(original[,sddistance], x[,"sddistance"]))
      } else {
          sddistance.val <- unique(c(x[,"sddistance"], 0.4175))
      }
      if (varying=="alpha") {
          alpha.val <- unique(c(original[,alpha], x[,"alpha"]))
      } else {
          alpha.val <- unique(x[,"alpha"])
      }
 output <- original[discount %in% discount.val& parties %in%parties.val  & alpha %in% alpha.val & sddistance %in% sddistance.val]
       output <- as.data.frame(output)
       if (by!="") {
        plot.data <- rbind(plot.data, 
                                cbind(x=output[,varying],
                                y=output[,interest.mean],
                               lower=output[,interest.mean]-(1.96*output[,interest.se]),
                               upper=output[,interest.mean]+(1.96*output[,interest.se]),
                                z=output[,by],
                                original=1)
                           )
#       plot.data <- as.data.frame(cbind(x,y, y.se, z))
       } else {    
        plot.data <- rbind(plot.data, 
                                cbind(x=output[,varying],
                                y=output[,interest.mean],
                               lower=output[,interest.mean]-(1.96*output[,interest.se]),
                               upper=output[,interest.mean]+(1.96*output[,interest.se]),
                                original=1)
                           )
      }
     }
     p <- ggplot(plot.data)
     p <- p + theme_bw()
     p <- p + theme(panel.grid.major.x = element_blank())
     p <- p + theme(panel.grid.major.y = element_blank())
     p <- p + theme(panel.grid.minor.x = element_blank())
     p <- p + theme(panel.grid.minor.y = element_blank())
     p <- p + xlab(varying)
     p <- p + ylab(interest.mean)
     if (by!="") {
      if (color==FALSE) {
       p <- p+ facet_grid(~z)
       p <- p + geom_ribbon(aes(x=x,
                                ymin=lower,
                                ymax=upper),
                            fill="grey74",
                            data=plot.data[plot.data$original==0,])
       p <- p + geom_line(aes(x=x, y=y),
                            data=plot.data[plot.data$original==0,])
      } else {
       p <- p + geom_ribbon(aes(x=x,
                                ymin=lower,
                                ymax=upper),
                            fill="grey74",
                            data=plot.data[plot.data$original==0,])
       p <- p + geom_line(aes(x=x, y=y, color=as.factor(z)),
                            data=plot.data[plot.data$original==0,])
      }
     } else {
         p <- p + geom_ribbon(aes(x=x,
                                 ymin=lower,
                                 ymax=upper),
                             fill="grey74",
                            data=plot.data[plot.data$original==0,])
         p <- p + geom_line(aes(x=x, y=y),
                            data=plot.data[plot.data$original==0,])
     }
     if (is.null(original)==FALSE) {
        p <- p + geom_segment(aes(x = x,
                                  y = lower,
                                  xend = x,
                                  yend = upper),
                              size=.1, data=plot.data[plot.data$original==1,])
       p <- p + geom_point(aes(x = x, y = y),
                           data=plot.data[plot.data$original==1,])
       p <- p+ facet_grid(~z)
     }
     return(p)
 } else {
     return(output)
 }
}


####################
# PLOT REAL VALUES #
####################

## parties:       3    4    5    6  7  8
## aspiration: 0.25 0.50 0.75 0.90
## memory:      2  3  4  5  6  7  8  9 10 
## sd.distance: 0.00 0.25 0.50 1.00 1.50 2.00
## discount:    0.00 0.10 0.25 0.50 0.75 0.90 1.00
## alpha:       0.00 0.10 0.25 0.50 0.75 0.90 1.00 (1=policy-seeking)

## plot.vars(ensemble, "parties", 
##           parties=5, discount=.5, alpha=.5, sd.distance=1)
## plot.vars(ensemble, "sd.distance", by="alpha", parties=5, discount=.5)
## plot.vars(ensemble, "discount","alpha", parties=5, discount=.5)
## plot.vars(ensemble, "alpha", parties=5, discount=.5, sd.distance=1)
## plot.vars(ensemble, "parties","alpha",
##           parties=5, discount=.5, alpha=.5, sd.distance=1)
## plot.vars(ensemble, "parties","sd.distance",
##           parties=5, discount=.5, alpha=.5, sd.distance=1)
## plot.vars(ensemble, "alpha", "sd.distance")



###########################################
# POLYWOG RESULTS, PREDICT, AND PLOT THEM #
###########################################


## cl <- makeCluster(4)
## registerDoParallel(cl)

## ##############
## # Aggregator #
## ##############
## print("Aggregator")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==2&ensemble$round==76,]
## formula=as.formula(mean.eccentricity~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.ecc, file="A.ecc.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.rep, file="A.rep.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.gov.rep, file="A.gov.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.seats~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.gov.seats <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.gov.seats, file="A.gov.seats.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.num.gov.parties~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.num.gov.parties <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.num.gov.parties, file="A.num.gov.parties.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.cohesion~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.gov.cohesion <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.gov.cohesion, file="A.gov.cohesion.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.oversized~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.gov.oversized <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.gov.oversized, file="A.oversized.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.majority~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.gov.majority <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.gov.majority, file="A.majority.Rdata")

## ## if (length(unique(data$mean.caretaker))!=1) {
## ##  set.seed(83758) #from random.org
## ##  formula=as.formula(mean.caretaker~
## ##                        alpha
## ##                        +parties
## ##                        +discount
## ##                        +sddistance)
## ##  A.caretaker <- polywog.select(formula=formula, data=data, degrees.cv=1:3, boot=50)
## ##  save(A.caretaker, file="A.caretaker.Rdata")
## ## }

## set.seed(83758) #from random.org
## formula=as.formula(mean.ENgov~
##                       alpha
##                       +parties
##                       +discount
##                       +sddistance)
## A.ENgov <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=40)
## save(A.ENgov, file="A.ENgov.Rdata")



## ##################
## # Sat Governator #
## ##################
## print("SatGov")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==3&ensemble$round==76,]
## formula=as.formula(mean.eccentricity~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.ecc, file="SG.ecc.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.rep, file="SG.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.gov.rep, file="SG.gov.rep.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.seats~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.gov.seats <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.gov.seats, file="SG.gov.seats.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.num.gov.parties~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.num.gov.parties <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.num.gov.parties, file="SG.num.gov.parties.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.cohesion~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.gov.cohesion <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.gov.cohesion, file="SG.gov.cohesion.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.oversized~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.gov.oversized <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.gov.oversized, file="SG.oversized.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.majority~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.gov.majority <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.gov.majority, file="SG.majority.Rdata")

## ## if (length(unique(data$mean.caretaker))!=1) {
## ##  set.seed(83758) #from random.org
## ##  formula=as.formula(mean.caretaker~
## ##                        alpha
## ##                        +parties
## ##                        +discount
## ##                        +sddistance)
## ##  SG.caretaker <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## ##  save(SG.caretaker, file="SG.caretaker.Rdata")
## ## }

## set.seed(83758) #from random.org
## formula=as.formula(mean.ENgov~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.ENgov <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.ENgov, file="SG.ENgov.Rdata")



## ###########
## # Sticker #
## ###########
## print("Sticker")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==4&ensemble$round==76,]
## formula=as.formula(mean.eccentricity~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.ecc, file="S.ecc.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.rep, file="S.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.gov.rep, file="S.gov.rep.Rdata")



## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.seats~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.gov.seats <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.gov.seats, file="S.gov.seats.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.num.gov.parties~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.num.gov.parties <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.num.gov.parties, file="S.num.gov.parties.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.cohesion~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.gov.cohesion <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.gov.cohesion, file="S.gov.cohesion.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.oversized~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.gov.oversized <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.gov.oversized, file="S.oversized.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.majority~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.gov.majority <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.gov.majority, file="S.majority.Rdata")

## ## if (length(unique(data$mean.caretaker))!=1) {
## ##  set.seed(83758) #from random.org
## ##  formula=as.formula(mean.caretaker~
## ##                        alpha
## ##                        +parties
## ##                        +discount
## ##                        +sddistance)
## ##  S.caretaker <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## ##  save(S.caretaker, file="S.caretaker.Rdata")
## ## }

## set.seed(83758) #from random.org
## formula=as.formula(mean.ENgov~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.ENgov <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.ENgov, file="S.ENgov.Rdata")




## ## HUNTER #

## print("Hunter")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==7]
## formula=as.formula(mean.eccentricity~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.ecc, file="H.ecc.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.rep, file="H.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.representation2~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.gov.rep, file="H.gov.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.seats~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.gov.seats <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.gov.seats, file="H.gov.seats.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.num.gov.parties~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.num.gov.parties <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.num.gov.parties, file="H.num.gov.parties.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.cohesion~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.gov.cohesion <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.gov.cohesion, file="H.gov.cohesion.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.oversized~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.gov.oversized <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.gov.oversized, file="H.oversized.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.majority~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.gov.majority <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.gov.majority, file="H.majority.Rdata")

## ## if (length(unique(data$mean.caretaker))!=1) {
## ##  set.seed(83758) #from random.org
## ##  formula=as.formula(mean.caretaker~
## ##                         alpha
## ##                        +parties
## ##                        +discount
## ##                        +sddistance)
## ##  H.caretaker <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## ##  save(H.caretaker, file="H.caretaker.Rdata")
## ## }

## set.seed(83758) #from random.org
## formula=as.formula(mean.ENgov~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.ENgov <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.ENgov, file="H.ENgov.Rdata")



## ###################
## # Hunt Governator #
## ###################
## print("HunGov")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==8&ensemble$round==76,]
## formula=as.formula(mean.eccentricity~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.ecc, file="HG.ecc.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.rep, file="HG.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.gov.rep, file="HG.gov.rep.Rdata")



## set.seed(83758) #from random.org
## formula=as.formula(mean.gov.seats~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.gov.seats <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.gov.seats, file="HG.gov.seats.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.num.gov.parties~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.num.gov.parties <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.num.gov.parties, file="HG.num.gov.parties.Rdata")


## ## gov.cohesion is always 0


## set.seed(83758) #from random.org
## formula=as.formula(mean.oversized~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.gov.oversized <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.gov.oversized, file="HG.oversized.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(mean.majority~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.gov.majority <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.gov.majority, file="HG.majority.Rdata")

## ## if (length(unique(data$mean.caretaker))!=1) {
## ##  set.seed(83758) #from random.org
## ##  formula=as.formula(mean.caretaker~
## ##                        alpha
## ##                        +parties
## ##                        +discount
## ##                        +sddistance)
## ##  HG.caretaker <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## ##  save(HG.caretaker, file="HG.caretaker.Rdata")
## ## }

## set.seed(83758) #from random.org
## formula=as.formula(mean.ENgov~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.ENgov <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.ENgov, file="HG.ENgov.Rdata")
## stopCluster(cl)



load("A.ecc.Rdata")
load("A.rep.Rdata")
load("A.gov.rep2.Rdata")
A.gov.rep2 <- A.gov.rep
load("A.gov.rep.Rdata")
load("A.num.gov.parties.Rdata")
load("A.gov.seats.Rdata")
load("A.gov.cohesion.Rdata")
load("A.majority.Rdata")
load("A.oversized.Rdata")
load("A.ENgov.Rdata")
load("SG.ecc.Rdata")
load("SG.rep.Rdata")
load("SG.gov.rep2.Rdata")
SG.gov.rep2 <- SG.gov.rep
load("SG.gov.rep.Rdata")
load("SG.num.gov.parties.Rdata")
load("SG.gov.seats.Rdata")
load("SG.gov.cohesion.Rdata")
load("SG.majority.Rdata")
load("SG.oversized.Rdata")
load("SG.ENgov.Rdata")
load("S.ecc.Rdata")
load("S.rep.Rdata")
load("S.gov.rep2.Rdata")
S.gov.rep2 <- S.gov.rep
load("S.gov.rep.Rdata")
load("S.num.gov.parties.Rdata")
load("S.gov.seats.Rdata")
load("S.gov.cohesion.Rdata")
load("S.majority.Rdata")
load("S.oversized.Rdata")
load("S.ENgov.Rdata")
load("H.ecc.Rdata")
load("H.rep.Rdata")
load("H.gov.rep2.Rdata")
H.gov.rep2 <- H.gov.rep
load("H.gov.rep.Rdata")
load("H.num.gov.parties.Rdata")
load("H.gov.seats.Rdata")
load("H.gov.cohesion.Rdata")
load("H.majority.Rdata")
load("H.oversized.Rdata")
load("H.ENgov.Rdata")
load("HG.ecc.Rdata")
load("HG.rep.Rdata")
load("HG.gov.rep2.Rdata")
HG.gov.rep2 <- HG.gov.rep
load("HG.gov.rep.Rdata")
load("HG.num.gov.parties.Rdata")
load("HG.gov.seats.Rdata")
## load("HG.gov.cohesion.Rdata") # is always 0
load("HG.majority.Rdata")
load("HG.oversized.Rdata")
load("HG.ENgov.Rdata")



v.parties <- seq(3,8,1)
v.discount <- .5
v.sddistance <- 1
v.alpha <- .5
Ae <- predict.polywog(A.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Se <- predict.polywog(S.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGe <- predict.polywog(HG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGe <- predict.polywog(SG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
He <- predict.polywog(H.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ar <- predict.polywog(A.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sr <- predict.polywog(S.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGr <- predict.polywog(HG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGr <- predict.polywog(SG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hr <- predict.polywog(H.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ag <- predict.polywog(A.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sg <- predict.polywog(S.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
SGg <- predict.polywog(SG.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hg <- predict.polywog(H.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
HGg <- predict.polywog(HG.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
t.data <- data.frame(x=c(7.5, 7.5, 7.6, 6.9, 4),
                     y=c(4.35, 3.85, 2.10, 1.17, .95),
                     label=c("S", "A", "H",
                         "SG", "Governator"),
                     type="eccentricity.mean")
data <- rbindlist(list(Ae, Se, He, SGe, HGe, Ar, Sr, HGr, SGr, Hr, Ag, Sg, Hg, SGg, HGg))
data$type <- as.character(data$type)
data$type <- ifelse(data$type=="mean.eccentricity", "Eccentricity", data$type)
data$type <- ifelse(data$type=="mean.representation", "System Misery", data$type)
data$type <- ifelse(data$type=="mean.gov.representation", "Government Misery", data$type)
data$type <- ifelse(data$type=="mean.gov.representation2", "Government Misery", data$type)
data$type <- factor(data$type, levels=c("Eccentricity", "System Misery", "Government Misery"))
data$means <- ifelse(data$type=="System Misery", data$means*(-1), data$means)
data$upper <- ifelse(data$type=="System Misery", data$upper*(-1), data$upper)
data$lower <- ifelse(data$type=="System Misery", data$lower*(-1), data$lower)
data <- data[,label.x:=max(parties)+.2, list(rule, type)]
data <- data[,label.y:=means[parties+.2==label.x], list(rule, type)]
data$label.y <- ifelse(data$type=="Eccentricity"&data$rule=="G", data$label.y+.001,data$label.y)
data$label.y <- ifelse(data$type=="System Misery"&data$rule=="S", data$label.y+.001,data$label.y)
data$label.y <- ifelse(data$type=="Government Misery"&data$rule=="S", data$label.y+.01,data$label.y)
p <- ggplot(data=data, aes(x=parties, y=means, group=rule))
p <- p + geom_ribbon(aes(ymin=lower, ymax=upper, x=parties), fill="gray80")
p <- p + facet_grid(type~., scales="free")
p <- p + geom_line(aes(colour=as.factor(rule)))
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank())
p <- p + theme(panel.grid.major.y = element_blank())
p <- p + theme(panel.grid.minor.x = element_blank())
p <- p + theme(panel.grid.minor.y = element_blank())
p <- p + ylab("Expected Mean of Mean")
p <- p + xlab("Number of Parties")
p <- p + geom_text(aes(x=label.x, y=label.y, label=rule))
## p <- p + scale_colour_manual(name="Rule", values=c("red", "gold2", "green", "blue", "magenta"))
p <- p + coord_fixed()
## p <- p + geom_dl(aes(label=as.factor(rule)), method = list(dl.trans(x = x + 0.2), "last.points", cex = 0.8))
p <- p + theme(legend.position="none")
#p <- p + scale_x_discrete(expand=c(0, .2))
p <- p + theme(text = element_text(size=15))
p
ggsave("parties.png", dpi=1000)



v.parties <- 5
v.discount <- .5
v.sddistance <- seq(0, 2, .1)
v.alpha <- c(0, .5, 1)
Ae <- predict.polywog(A.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Se <- predict.polywog(S.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGe <- predict.polywog(HG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGe <- predict.polywog(SG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
He <- predict.polywog(H.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ar <- predict.polywog(A.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sr <- predict.polywog(S.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGr <- predict.polywog(HG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGr <- predict.polywog(SG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hr <- predict.polywog(H.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ag <- predict.polywog(A.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sg <- predict.polywog(S.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
SGg <- predict.polywog(SG.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hg <- predict.polywog(H.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
HGg <- predict.polywog(HG.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
data <- rbindlist(list(Ag, Sg, Hg, SGg, HGg))
data$alpha <- factor(data$alpha*2, labels=c("100% office-motivated",
                                    "50% office-motivated\n50% policy-motivated",
                                    "100% policy-motivated"))
data <- data[,label.x:=max(sddistance)+.2, list(rule, alpha)]
data <- data[,label.y:=means[sddistance+.2==label.x], list(rule, alpha)]
data$label.y <- ifelse(data$alpha!="100% policy-motivated"&data$rule=="G",
                       data$label.y+.015,data$label.y)
p <- ggplot(data=data, aes(x=sddistance, y=means, group=rule))
p <- p + geom_ribbon(aes(ymin=lower, ymax=upper, x=sddistance), fill="gray80")
p <- p + facet_grid(~alpha, scales="free")
p <- p + geom_line(aes(colour=as.factor(rule)))
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank())
p <- p + theme(panel.grid.major.y = element_blank())
p <- p + theme(panel.grid.minor.x = element_blank())
p <- p + theme(panel.grid.minor.y = element_blank())
p <- p + theme(legend.position="bottom")
p <- p + ylab("Expected Mean of Government Representativeness")
p <- p + xlab("Dispersion of Parties' Ideal Positions")
p <- p + scale_colour_discrete(name="Rule")
p <- p + coord_fixed()
p <- p + theme(legend.position="none")
p <- p + scale_x_discrete(expand=c(0, .2))
p <- p + theme(text = element_text(size=15))
p <- p + geom_text(aes(x=label.x, y=label.y, label=rule))
p
ggsave("sddistance.png", dpi=1000)



v.parties <- 5
v.discount <- seq(0, 1, .01)
v.sddistance <- 1
v.alpha <- c(0, .5, 1)
Ae <- predict.polywog(A.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Se <- predict.polywog(S.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGe <- predict.polywog(HG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGe <- predict.polywog(SG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
He <- predict.polywog(H.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ar <- predict.polywog(A.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sr <- predict.polywog(S.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGr <- predict.polywog(HG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGr <- predict.polywog(SG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hr <- predict.polywog(H.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ag <- predict.polywog(A.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sg <- predict.polywog(S.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
SGg <- predict.polywog(SG.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hg <- predict.polywog(H.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
HGg <- predict.polywog(HG.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
t.data <- data.frame(x=c(7.5, 7.5, 7.6, 6.9, 4),
                     y=c(4.35, 3.85, 2.10, 1.17, .95),
                     label=c("S", "A", "H",
                         "SG", "G"),
                     type="eccentricity.mean")
data <- rbindlist(list(Ag, Sg, Hg, SGg, HGg))
data$alpha <- factor(data$alpha*2, labels=c("100% office-motivated",
                                    "50% office-motivated\n50% policy-motivated",
                                    "100% policy-motivated"))
data <- data[,label.x:=max(sddistance)+.2, list(rule, alpha)]
data <- data[,label.y:=means[sddistance+.2==label.x], list(rule, alpha)]
## data$label.y <- ifelse(data$alpha!="100% policy-motivated"&data$rule=="G",
##                        data$label.y+.015,data$label.y)
p <- ggplot(data=data, aes(x=discount, y=means, group=rule))
p <- p + geom_ribbon(aes(ymin=lower, ymax=upper, x=discount), fill="gray80")
p <- p + facet_grid(~alpha, scales="free")
p <- p + geom_line(aes(colour=as.factor(rule)))
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank())
p <- p + theme(panel.grid.major.y = element_blank())
p <- p + theme(panel.grid.minor.x = element_blank())
p <- p + theme(panel.grid.minor.y = element_blank())
p <- p + ylab("Expected Mean of Government Representativeness")
p <- p + xlab("Discount Factor for Caretaker Governments")
p <- p + scale_colour_discrete(name="Rule")
p <- p + coord_fixed()
p <- p + theme(legend.position="none")
p <- p + scale_x_discrete(expand=c(0, .2))
p <- p + theme(text = element_text(size=15))
p <- p + geom_text(aes(x=label.x, y=label.y, label=rule))
p
ggsave("discount.png", dpi=2000)

v.parties <- 5
v.discount <- .5
v.sddistance <- 1
v.alpha <- seq(0, 1, .01)
Ae <- predict.polywog(A.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Se <- predict.polywog(S.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGe <- predict.polywog(HG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGe <- predict.polywog(SG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
He <- predict.polywog(H.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ar <- predict.polywog(A.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sr <- predict.polywog(S.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGr <- predict.polywog(HG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGr <- predict.polywog(SG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hr <- predict.polywog(H.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ag <- predict.polywog(A.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sg <- predict.polywog(S.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
SGg <- predict.polywog(SG.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hg <- predict.polywog(H.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
HGg <- predict.polywog(HG.gov.rep2, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
t.data <- data.frame(x=c(7.5, 7.5, 7.6, 6.9, 4),
                     y=c(4.35, 3.85, 2.10, 1.17, .95),
                     label=c("S", "A", "H",
                         "SG", "G"),
                     type="eccentricity.mean")
data <- rbindlist(list(Ag, Sg, Hg, SGg, HGg))
data <- data[,label.x:=max(alpha)+.03, list(rule)]
data <- data[,label.y:=means[alpha+.03==label.x], list(rule)]
data$label.y <- ifelse(data$rule=="G",
                       data$label.y+.005,data$label.y)
p <- ggplot(data=data, aes(x=alpha, y=means, group=rule))
p <- p + geom_ribbon(aes(ymin=lower, ymax=upper, x=alpha), fill="gray80")
p <- p + geom_line(aes(colour=as.factor(rule)))
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank())
p <- p + theme(panel.grid.major.y = element_blank())
p <- p + theme(panel.grid.minor.x = element_blank())
p <- p + theme(panel.grid.minor.y = element_blank())
p <- p + ylab("Expected Mean of Government Representativeness")
p <- p + xlab("Level of Policy-Seekingness")
p <- p + theme(legend.position="none")
p <- p + scale_x_discrete(expand=c(0, .03))
p <- p + theme(text = element_text(size=15))
p <- p + geom_text(aes(x=label.x, y=label.y, label=rule))
p
ggsave("alpha.png", dpi=1000)



##############################


#Make Predictions based on data
ecc <- predict.polywog(A.ecc, parties=seq(3,8,1), discount=.5, sddistance=1, alpha=seq(0,1,.25))
rep <- predict.polywog(A.rep, parties=seq(3,8,1), discount=.5, sddistance=1, alpha=seq(0,1,.25))
gov <- predict.polywog(A.gov.rep, parties=seq(3,8,1), discount=.5, sddistance=1, alpha=seq(0,1,.25))

#Plot based on predictions which are based on the data-mined model.
sub <- ensemble[rule==2]
plot.polywog(ecc, varying="parties", by="alpha", original=sub, interest="eccentricity")
plot.polywog(rep, varying="parties", by="alpha", original=sub, interest="representation")
plot.polywog(gov, varying="parties", by="alpha", original=sub, interest="rep.gov")



rep <- predict.polywog(SG.rep, parties=seq(3,8,1), discount=.5, sddistance=1, alpha=seq(0,1,.25))
sub <- ensemble[rule==3]
plot.polywog(rep, varying="parties", by="alpha", original=sub, interest="representation")

rep <- predict.polywog(SG.rep, parties=seq(3,8,1), discount=.5,
                       sddistance=unique(ensemble[,sddistance]), alpha=.5)
sub <- ensemble[rule==3]
plot.polywog(rep, varying="sddistance", by="parties", original=sub, interest="representation")

rep <- predict.polywog(SG.rep, parties=seq(3,8,1), discount=unique(ensemble[,discount]),
                       sddistance=1, alpha=.5)
sub <- ensemble[rule==3]
plot.polywog(rep, varying="discount", by="parties", original=sub, interest="representation")



rep <- predict.polywog(HG.rep, parties=seq(3,8,1), discount=.5, sddistance=1, alpha=seq(0,1,.25))
sub <- ensemble[rule==8]
plot.polywog(rep, varying="parties", by="alpha", original=sub, interest="representation")

rep <- predict.polywog(HG.rep, parties=seq(3,8,1), discount=.5,
                       sddistance=unique(ensemble[,sddistance]), alpha=.5)
sub <- ensemble[rule==8]
plot.polywog(rep, varying="sddistance", by="parties", original=sub, interest="representation")

rep <- predict.polywog(HG.rep, parties=seq(3,8,1), discount=unique(ensemble[,discount]),
                       sddistance=1, alpha=.5)
sub <- ensemble[rule==8]
plot.polywog(rep, varying="discount", by="parties", original=sub, interest="representation")


ecc <- predict.polywog(SG.ecc, parties=seq(3,8,1), discount=.5, sddistance=1, alpha=seq(0,1,.25))
sub <- ensemble[rule==3]
plot.polywog(ecc, varying="parties", by="alpha", original=sub, interest="eccentricity")

ecc <- predict.polywog(SG.ecc, parties=seq(3,8,1), discount=unique(ensemble[,discount]), sddistance=1, alpha=1)
sub <- ensemble[rule==3]
plot.polywog(ecc, varying="parties", by="discount", original=sub, interest="eccentricity")


######################################################################################

## v.parties <- unique(ensemble[,parties])
## v.alpha <- unique(ensemble[,alpha])
## v.discount <- .5
## v.sddistance <- 1
## varying <- "parties"
## by <- "alpha"
## original <- ensemble[round==76]

complete <- function(v.parties, v.alpha, v.discount, v.sddistance,
                     varying, by, original,
                     A.ecc, A.rep, A.gov.rep,
                     SG.ecc, SG.rep, SG.gov.rep,
                     S.ecc, S.rep, S.gov.rep,
                     H.ecc, H.rep, H.gov.rep,
                     HG.ecc, HG.rep, HG.gov.rep) {
Ae <- predict.polywog(A.ecc, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Aggregator")
SGe <- predict.polywog(SG.ecc, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Satisficing\nGovernator")
Se <- predict.polywog(S.ecc, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Sticker")
He <- predict.polywog(H.ecc, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Hunter")
HGe <- predict.polywog(HG.ecc, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Governator")
Ar <- predict.polywog(A.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Aggregator")
SGr <- predict.polywog(SG.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Satisficing\nGovernator")
Sr <- predict.polywog(S.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Sticker")
Hr <- predict.polywog(H.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Hunter")
HGr <- predict.polywog(HG.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Governator")
Ag <- predict.polywog(A.gov.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Aggregator")
SGg <- predict.polywog(SG.gov.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Satisficing\nGovernator")
Sg <- predict.polywog(S.gov.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Sticker")
Hg <- predict.polywog(H.gov.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Hunter")
HGg <- predict.polywog(HG.gov.rep, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Governator")
data <- rbindlist(list(Ae, Se, He, SGe, HGe,
                       Ar, Sr, Hr, SGr, HGr,
                       Ag, Sg, Hg, SGg, HGg))
data$mean.eccentricity <- NULL
data$mean.representation <- NULL
data$mean.gov.representation <- NULL
if (is.data.frame(original)==TRUE) {
    d <- original[discount %in% v.discount
                 & parties %in% v.parties
                 & alpha %in% v.alpha
                 & sddistance %in% v.sddistance,
                  list(rule, alpha, parties, discount, sddistance,
                    mean.eccentricity, se.eccentricity,
                    mean.representation, se.representation,
                    mean.gov.representation, se.gov.representation)]
  d$rule <- ifelse(d$rule==2, "Aggregator",
                   ifelse(d$rule==3, "Satisficing\nGovernator",
                          ifelse(d$rule==4, "Sticker",
                                 ifelse(d$rule==7, "Hunter", "Governator"))))
  d$rule2 <- factor(d$rule, labels=levels(data$rule))
  setnames(d, c("mean.eccentricity", "se.eccentricity",
                "mean.representation", "se.representation",
                "mean.gov.representation", "se.gov.representation"),
           c("mean1", "se1", "mean2", "se2", "mean3", "se3"))
  d <- melt(d, measure=patterns("mean.", "se."), value.name=c("Tmeans", "Tse"),
            variable.name="type")
  d$type <- factor(d$type,
                   labels=c("mean.eccentricity", "mean.representation", "mean.gov.representation"))
  d$Tlower <- d$Tmeans-(1.96*d$Tse)
  d$Tupper <- d$Tmeans+(1.96*d$Tse)
  d$Tse <- NULL
  setkey(d, rule, alpha, parties, discount, sddistance, type)
  setkey(data, rule, alpha, parties, discount, sddistance, type)
  data <- data[d]    
}
 data <- as.data.frame(data)
 plot.data <-data.frame(x=data[,varying],
                          y=data$means,
                          lower=data$lower,
                          upper=data$upper,
                          rule=data$rule,
                          type=data$type)
 if (by!="") {
   plot.data$z <- data[,by]
 }
 if (is.data.frame(original)==TRUE) {
   plot.data$Tmeans <- data$Tmeans
   plot.data$Tlower <- data$Tlower
   plot.data$Tupper <- data$Tupper
 }
     p <- ggplot(plot.data, aes(x=x, y=y, group=rule))
     #p <- p + geom_segment(aes(colour=variable, y=upper, yend=lower, xend=x))
     #p <- p + geom_point(aes(colour=variable))
     p <- p + geom_ribbon(fill="grey73", aes(ymin=lower, ymax=upper))
     p <- p + geom_line(aes(colour=as.factor(rule)))
     p <- p + theme_bw()
     p <- p + theme(panel.grid.major.x = element_blank())
     p <- p + theme(panel.grid.major.y = element_blank())
     p <- p + theme(panel.grid.minor.x = element_blank())
     p <- p + theme(panel.grid.minor.y = element_blank())
     p <- p + ylab("")
     p <- p + scale_colour_manual("Rules",
                                  values=c("red", "deeppink", "green", "blue", "gold"),
                                  labels=c("Aggregator", "Governator", "Hunter",
                                        "Satisficing\nGovernator", "Sticker"))
     if (is.data.frame(original)==TRUE) {
      p <- p + geom_point(aes(x=x, y=Tmeans, colour=as.factor(rule)))
      p <- p + geom_segment(aes(xend=x, y=Tlower, yend=Tupper, colour=as.factor(rule)))
     }
     if (by!="") {
         p <- p + facet_grid(type~z, scales="free")
         p <- p + xlab(paste0(varying, " by ", by))
      } else {
         p <- p + facet_grid(type~., scales="free")
         p <- p + xlab(varying)
      }
return(p)   
}

pdf("gov.rep.pdf")
###########
# parties #
###########
complete(v.parties=3:8, v.alpha=unique(ensemble[,alpha]),
         v.discount=.5, v.sddistance=1,
         varying="parties", by="alpha", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

complete(v.parties=3:8, v.alpha=.5,
         v.discount=unique(ensemble[,discount]), v.sddistance=1,
         varying="parties", by="discount", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

complete(v.parties=3:8, v.alpha=.5,
         v.discount=.5, v.sddistance=unique(ensemble[,sddistance]),
         varying="parties", by="sddistance", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

#########
# alpha #
#########
complete(v.parties=3:8, v.alpha=unique(ensemble[,alpha]),
         v.discount=.5, v.sddistance=1,
         varying="alpha", by="parties", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

complete(v.parties=4, v.alpha=unique(ensemble[,alpha]),
         v.discount=unique(ensemble[,discount]), v.sddistance=1,
         varying="alpha", by="discount", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

complete(v.parties=4, v.alpha=unique(ensemble[,alpha]),
         v.discount=.5, v.sddistance=unique(ensemble[,sddistance]),
         varying="alpha", by="sddistance", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

###########
#discount #
###########
complete(v.parties=3:8, v.alpha=.5,
         v.discount=unique(ensemble[,discount]), v.sddistance=1,
         varying="discount", by="parties", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

complete(v.parties=4, v.alpha=unique(ensemble[,alpha]),
         v.discount=unique(ensemble[,discount]), v.sddistance=1,
         varying="discount", by="alpha", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

complete(v.parties=4, v.alpha=.5,
         v.discount=unique(ensemble[,discount]), v.sddistance=unique(ensemble[,sddistance]),
         varying="discount", by="sddistance", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)


##############
# sddistance #
##############
complete(v.parties=3:8, v.alpha=.5,
         v.discount=.5, v.sddistance=unique(ensemble[,discount]),
         varying="sddistance", by="parties", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

complete(v.parties=4, v.alpha=unique(ensemble[,alpha]),
         v.discount=.5, v.sddistance=unique(ensemble[,discount]),
         varying="sddistance", by="alpha", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)

complete(v.parties=4, v.alpha=.5,
         v.discount=unique(ensemble[,discount]), v.sddistance=unique(ensemble[,discount]),
         varying="sddistance", by="discount", original=ensemble[round==76],
         A.ecc=A.ecc, A.rep=A.rep, A.gov.rep=A.gov.rep,
         S.ecc=S.ecc, S.rep=S.rep, S.gov.rep=S.gov.rep,
         H.ecc=H.ecc, H.rep=H.rep, H.gov.rep=H.gov.rep,
         HG.ecc=HG.ecc, HG.rep=HG.rep, HG.gov.rep=HG.gov.rep,
         SG.ecc=SG.ecc, SG.rep=SG.rep, SG.gov.rep=SG.gov.rep)
dev.off()

##############################################################################################



## v.parties <- unique(ensemble[,parties])
## v.alpha <- unique(ensemble[,alpha])
## v.discount <- .5
## v.sddistance <- 1
## varying <- "parties"
## by <- "alpha"
## original <- ensemble[round==76]

complete.gov <- function(v.parties, v.alpha, v.discount, v.sddistance,
                         varying, by, original,
 A.gov.seats, A.gov.cohesion, A.num.gov.parties, A.oversized, A.majority, A.ENgov,
 SG.gov.seats, SG.gov.cohesion, SG.num.gov.parties, SG.oversized,SG.majority, SG.ENgov,
 S.gov.seats, S.gov.cohesion, S.num.gov.parties, S.oversized, S.majority, S.ENgov,
 H.gov.seats, H.gov.cohesion, H.num.gov.parties, H.oversized, H.majority, H.ENgov,
 HG.gov.seats, HG.gov.cohesion, HG.num.gov.parties, HG.oversized, HG.majority, HG.ENgov) {
Aseats <- predict.polywog(A.gov.seats, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Aggregator")
Anum <- predict.polywog(A.num.gov.parties, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Aggregator")
Acohesion <- predict.polywog(A.gov.cohesion, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Aggregator")
Aoversized <- predict.polywog(A.gov.oversized, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Aggregator")
Amajority <- predict.polywog(A.gov.majority, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Aggregator")
AENgov <- predict.polywog(A.ENgov, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Aggregator")
SGseats <- predict.polywog(SG.gov.seats, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Satisficing\nGovernator")
SGnum <- predict.polywog(SG.num.gov.parties, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Satisficing\nGovernator")
SGcohesion <- predict.polywog(SG.gov.cohesion, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Satisficing\nGovernator")
SGoversized <- predict.polywog(SG.gov.oversized, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Satisficing\nGovernator")
SGmajority <- predict.polywog(SG.gov.majority, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Satisficing\nGovernator")
SGENgov <- predict.polywog(SG.ENgov, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Satisficing\nGovernator")
HGseats <- predict.polywog(HG.gov.seats, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Governator")
HGnum <- predict.polywog(HG.num.gov.parties, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Governator")
## HGcohesion <- predict.polywog(HG.gov.cohesion, parties=v.parties, discount=v.discount,
##                       sddistance=v.sddistance, alpha=v.alpha, rule="Governator")
HGcohesion <- as.data.table(expand.grid(parties=v.parties, alpha=v.alpha, discount=v.discount,
                            sddistance=v.sddistance, lower=0, upper=0, means=0,
                         type="mean.gov.cohesion", rule="Governator"))
HGoversized <- predict.polywog(HG.gov.oversized, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Governator")
HGmajority <- predict.polywog(HG.gov.majority, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Governator")
HGENgov <- predict.polywog(HG.ENgov, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Governator")
Sseats <- predict.polywog(S.gov.seats, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Sticker")
Snum <- predict.polywog(S.num.gov.parties, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Sticker")
Scohesion <- predict.polywog(S.gov.cohesion, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Sticker")
Soversized <- predict.polywog(S.gov.oversized, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Sticker")
Smajority <- predict.polywog(S.gov.majority, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Sticker")
SENgov <- predict.polywog(S.ENgov, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Sticker")
Hseats <- predict.polywog(H.gov.seats, parties=v.parties, discount=v.discount,
                     sddistance=v.sddistance, alpha=v.alpha, rule="Hunter")
Hnum <- predict.polywog(H.num.gov.parties, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Hunter")
Hcohesion <- predict.polywog(H.gov.cohesion, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Hunter")
Hoversized <- predict.polywog(H.gov.oversized, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Hunter")
Hmajority <- predict.polywog(H.gov.majority, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Hunter")
HENgov <- predict.polywog(H.ENgov, parties=v.parties, discount=v.discount,
                      sddistance=v.sddistance, alpha=v.alpha, rule="Hunter")
data <- rbindlist(list(Aseats, Sseats, Hseats, SGseats, HGseats,
                       Anum, Snum, Hnum, SGnum, HGnum,
                       Acohesion, Scohesion, Hcohesion, SGcohesion, HGcohesion,
                       Aoversized, Soversized, Hoversized, SGoversized, HGoversized,
                       Amajority, Smajority, Hmajority, SGmajority, HGmajority,
                       AENgov, SENgov, HENgov, SGENgov, HGENgov), fill=TRUE)
data$mean.eccentricity <- NULL
data$mean.representation <- NULL
data$mean.gov.representation <- NULL
data$mean.gov.seats <- NULL
data$mean.num.gov.parties <- NULL
data$mean.gov.cohesion <- NULL
data$mean.oversized <- NULL
data$mean.majority <- NULL
data$mean.ENgov <- NULL
if (is.data.frame(original)==TRUE) {
    d <- original[discount %in% v.discount
                 & parties %in% v.parties
                 & alpha %in% v.alpha
                 & sddistance %in% v.sddistance,
                  list(rule, alpha, parties, discount, sddistance,
                    mean.gov.seats, se.gov.seats,
                    mean.num.gov.parties, se.num.gov.parties,
                    mean.gov.cohesion, se.gov.cohesion,
                    mean.oversized, se.oversized,
                    mean.majority, se.majority,
                    mean.ENgov, se.ENgov)]
  d$rule <- ifelse(d$rule==2, "Aggregator",
                   ifelse(d$rule==3, "Satisficing\nGovernator",
                          ifelse(d$rule==4, "Sticker",
                                 ifelse(d$rule==7, "Hunter", "Governator"))))
  d$rule2 <- factor(d$rule, labels=levels(data$rule))
  setnames(d, c(    "mean.gov.seats", "se.gov.seats",
                    "mean.num.gov.parties", "se.num.gov.parties",
                    "mean.gov.cohesion", "se.gov.cohesion",
                    "mean.oversized", "se.oversized",
                    "mean.majority", "se.majority",
                    "mean.ENgov", "se.ENgov"),
           c("mean1", "se1", "mean2", "se2", "mean3", "se3", "mean4", "se4", "mean5", "se5",
             "mean6", "se6"))
  d <- melt(d, measure=patterns("mean.", "se."), value.name=c("Tmeans", "Tse"),
            variable.name="type")
  d$type <- factor(d$type,
                   labels=c("mean.gov.seats", "mean.num.gov.parties", "mean.gov.cohesion",
                             "mean.oversized", "mean.majority", "mean.ENgov"))
  d$Tlower <- d$Tmeans-(1.96*d$Tse)
  d$Tupper <- d$Tmeans+(1.96*d$Tse)
  d$Tse <- NULL
  setkey(d, rule, alpha, parties, discount, sddistance, type)
  setkey(data, rule, alpha, parties, discount, sddistance, type)
  data <- data[d]    
}
 data <- as.data.frame(data)
 plot.data <-data.frame(x=data[,varying],
                          y=data$means,
                          lower=data$lower,
                          upper=data$upper,
                          rule=data$rule,
                          type=data$type)
 if (by!="") {
   plot.data$z <- data[,by]
 }
 if (is.data.frame(original)==TRUE) {
   plot.data$Tmeans <- data$Tmeans
   plot.data$Tlower <- data$Tlower
   plot.data$Tupper <- data$Tupper
 }
     p <- ggplot(plot.data, aes(x=x, y=y, group=rule))
     #p <- p + geom_segment(aes(colour=variable, y=upper, yend=lower, xend=x))
     #p <- p + geom_point(aes(colour=variable))
     p <- p + geom_ribbon(fill="grey73", aes(ymin=lower, ymax=upper))
     p <- p + geom_line(aes(colour=as.factor(rule)))
     p <- p + theme_bw()
     p <- p + theme(panel.grid.major.x = element_blank())
     p <- p + theme(panel.grid.major.y = element_blank())
     p <- p + theme(panel.grid.minor.x = element_blank())
     p <- p + theme(panel.grid.minor.y = element_blank())
     p <- p + ylab("")
     p <- p + scale_colour_manual("Rules",
                                  values=c("red", "deeppink", "green", "blue", "gold"),
                                  labels=c("Aggregator", "Governator", "Hunter",
                                        "Satisficing\nGovernator", "Sticker"))
     if (is.data.frame(original)==TRUE) {
      p <- p + geom_point(aes(x=x, y=Tmeans, colour=as.factor(rule)))
      p <- p + geom_segment(aes(xend=x, y=Tlower, yend=Tupper, colour=as.factor(rule)))
     }
     if (by!="") {
         p <- p + facet_grid(type~z, scales="free")
         p <- p + xlab(paste0(varying, " by ", by))
      } else {
         p <- p + facet_grid(type~., scales="free")
         p <- p + xlab(varying)
      }
return(p)   
}

pdf("governmentts.pdf")
###########
# parties #
###########
complete.gov(v.parties=3:8, v.alpha=unique(ensemble[,alpha]),
         v.discount=.5, v.sddistance=1,
         varying="parties", by="alpha", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)
             

complete.gov(v.parties=3:8, v.alpha=.5,
         v.discount=unique(ensemble[,discount]), v.sddistance=1,
         varying="parties", by="discount", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)

complete.gov(v.parties=3:8, v.alpha=.5,
         v.discount=.5, v.sddistance=unique(ensemble[,sddistance]),
         varying="parties", by="sddistance", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)

#########
# alpha #
#########
complete.gov(v.parties=3:8, v.alpha=unique(ensemble[,alpha]),
         v.discount=.5, v.sddistance=1,
         varying="alpha", by="parties", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)

complete.gov(v.parties=4, v.alpha=unique(ensemble[,alpha]),
         v.discount=unique(ensemble[,discount]), v.sddistance=1,
         varying="alpha", by="discount", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)

complete.gov(v.parties=4, v.alpha=unique(ensemble[,alpha]),
         v.discount=.5, v.sddistance=unique(ensemble[,sddistance]),
         varying="alpha", by="sddistance", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)

###########
#discount #
###########
complete.gov(v.parties=3:8, v.alpha=.5,
         v.discount=unique(ensemble[,discount]), v.sddistance=1,
         varying="discount", by="parties", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)

complete.gov(v.parties=4, v.alpha=unique(ensemble[,alpha]),
         v.discount=unique(ensemble[,discount]), v.sddistance=1,
         varying="discount", by="alpha", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)

complete.gov(v.parties=4, v.alpha=.5,
         v.discount=unique(ensemble[,discount]), v.sddistance=unique(ensemble[,sddistance]),
         varying="discount", by="sddistance", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)


##############
# sddistance #
##############
complete.gov(v.parties=3:8, v.alpha=.5,
         v.discount=.5, v.sddistance=unique(ensemble[,discount]),
         varying="sddistance", by="parties", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)

complete.gov(v.parties=4, v.alpha=unique(ensemble[,alpha]),
         v.discount=.5, v.sddistance=unique(ensemble[,discount]),
         varying="sddistance", by="alpha", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)

complete.gov(v.parties=4, v.alpha=.5,
         v.discount=unique(ensemble[,discount]), v.sddistance=unique(ensemble[,discount]),
         varying="sddistance", by="discount", original=ensemble[round==76],
          A.gov.seats=A.gov.seats, A.gov.cohesion=A.gov.cohesion,
          A.num.gov.parties=A.num.gov.parties, A.oversize=A.oversized, A.majority=A.majority,
          A.ENgov=A.ENgov,
          S.gov.seats=S.gov.seats, S.gov.cohesion=S.gov.cohesion,
          S.num.gov.parties=S.num.gov.parties, S.oversize=S.oversized, S.majority=S.majority,
          S.ENgov=S.ENgov,
          HG.gov.seats=HG.gov.seats, HG.gov.cohesion=HG.gov.cohesion,
          HG.num.gov.parties=HG.num.gov.parties, HG.oversize=HG.oversized, HG.majority=HG.majority,
          HG.ENgov=HG.ENgov,
          H.gov.seats=H.gov.seats, H.gov.cohesion=H.gov.cohesion,
          H.num.gov.parties=H.num.gov.parties, H.oversize=H.oversized, H.majority=H.majority,
          H.ENgov=H.ENgov,
          SG.gov.seats=SG.gov.seats, SG.gov.cohesion=SG.gov.cohesion,
          SG.num.gov.parties=SG.num.gov.parties, SG.oversize=SG.oversized, SG.majority=SG.majority,
          SG.ENgov=SG.ENgov)
dev.off()

##################################################################################################

v.parties <- seq(3,8,1)
v.discount <- .5
v.sddistance <- 1
v.alpha <- c(0, .5, 1)
Acohesion <- predict.polywog(A.gov.cohesion,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Aggregator")
Scohesion <- predict.polywog(S.gov.cohesion,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Sticker")
Hcohesion <- predict.polywog(H.gov.cohesion,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Hunter")
SGcohesion <- predict.polywog(SG.gov.cohesion,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Satsificing\nGovernator")
HGcohesion <- as.data.table(expand.grid(parties=v.parties, alpha=v.alpha, discount=v.discount,
                            sddistance=v.sddistance, lower=0, upper=0, means=0,
                         type="mean.gov.cohesion", rule="Governator"))
Hcohesion <- predict.polywog(H.gov.cohesion,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Hunter")
Amajority <- predict.polywog(A.gov.majority,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Aggregator")
Smajority <- predict.polywog(S.gov.majority,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Sticker")
HGmajority <- predict.polywog(HG.gov.majority,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Governator")
SGmajority <- predict.polywog(SG.gov.majority,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Satsificing\nGovernator")
Hmajority <- predict.polywog(H.gov.majority,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Hunter")
AENgov <- predict.polywog(A.ENgov,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Aggregator")
SENgov <- predict.polywog(S.ENgov,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Sticker")
HGENgov <- predict.polywog(HG.ENgov,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Governator")
SGENgov <- predict.polywog(SG.ENgov,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Satsificing\nGovernator")
HENgov <- predict.polywog(H.ENgov,
                             parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                             alpha=v.alpha, rule="Hunter")
data <- rbindlist(list(Acohesion, Scohesion, HGcohesion, SGcohesion, Hcohesion,
                       Amajority, Smajority, HGmajority, SGmajority, Hmajority,
                       AENgov, SENgov, HGENgov, SGENgov, HENgov), fill=TRUE)
data$mean.eccentricity <- NULL
data$mean.representation <- NULL
data$mean.gov.representation <- NULL
data$mean.gov.seats <- NULL
data$mean.num.gov.parties <- NULL
data$mean.gov.cohesion <- NULL
data$mean.oversized <- NULL
data$mean.majority <- NULL
data$mean.ENgov <- NULL
levels(data$type) <- c("Government\nCohesion", "Majority\nGovernment", "Effective Num\nof Parties")
data$alpha <- as.factor(data$alpha)
levels(data$alpha) <- c("100%\noffice-motivated", "50%office-motivated\n50%policy-motivated",
                        "100%\npolicy-motivated")
p <- ggplot(data=data, aes(x=parties, y=means, group=rule))
p <- p + geom_ribbon(aes(ymin=lower, ymax=upper, x=parties), fill="gray80")
p <- p + facet_grid(type~alpha, scales="free")
p <- p + geom_line(aes(colour=as.factor(rule)))
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank())
p <- p + theme(panel.grid.major.y = element_blank())
p <- p + theme(panel.grid.minor.x = element_blank())
p <- p + theme(panel.grid.minor.y = element_blank())
p <- p + ylab("Expected Mean of Mean")
p <- p + xlab("Number of Parties")
p <- p + scale_colour_manual(name="Rule", values=c("red", "gold2", "magenta", "blue", "green"))
p
ggsave("cohesion.png", dpi=200)
